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ABSTRACT 

We study the sample variance of the matter power spectrum for the standard A Cold Dark Matter 
universe. We use a total of 5000 cosmological TV-body simulations to study in detail the distribution of 
best-fit cosmological parameters and the baryon acoustic peak positions. The obtained distribution is 
compared with the results from the Fisher matrix analysis with and without including non-Gaussian 
errors. For the Fisher matrix analysis, we compute the derivatives of the matter power spectrum with 
respect to cosmological parameters using directly full nonlinear simulations. We show that the non- 
Gaussian errors increase the unmarginalized errors by up to a factor 5 for /cmax = 0.4/i/Mpc if there 
is only one free parameter provided other parameters are well determined by external information. 
On the other hand, for multi-parameter fitting, the impact of the non-Gaussian errors is significantly 
mitigated due to severe parameter degeneracies in the power spectrum. The distribution of the 
acoustic peak positions is well described by a Gaussian distribution, with its width being consistent 
with the statistical interval predicted from the Fisher matrix. We also examine systematic bias in 
the best-fit parameter due to the non-Gaussian errors. The bias is found to be smaller than the Icr 
statistical error for both the cosmological parameters and the acoustic scale positions. 
Subject headings: cosmology: theory - large-scale structure of universe 



1. INTRODUCTION 

The baryon acoustic oscillation (BAO) is imprinted in 
the distribution of galaxies as is found in the temperature 
fluctuations in the cosmic microwave background. The 
acoustic length scale is determined by the sound hori- 
zon of the photon-baryon fluid at recombination epoch; 
it can thus be used as a standard ruler which provides us 
with a robust method to measure distance scales out to 
essentially any epoch (e.g., Eisenstein et al. 1998; Blake 
& Glazebrook 2003; Seo & Eisenstein 2003; Matsubara 
2004; Guzik et al. 2007). Using the observed distance- 
redshift relation, we can obtain an accurate cosmic ex- 
pansion history, which in turn gives strong constraints on 
the nature of dark energy. The large-area galaxy surveys 
such as two-degree Field Survey (2dF) and Sloan Digital 
Sky Survey (SDSS) detected the BAO signature in the 
galaxy distribution (Cole et al. 2005; Eisenstein et al 
2005; Percival et al. 2007; Okumura et al. 2008; Gaz- 
tanagaetal. 2008; Sanchez et al. 2009). The latest result 
of the SDSS DR7 showed a constraint on the distance to 
a redshift z = 0.28 within 2.7% accuracy (Percival et al. 
2009; Reid et al. 2009; Kazin et al. 2009). Future and 
ongoing surveys such as the Baryon Oscillation Spectro- 
scopic Survey (BOSS)^, the Hobby-Eberly Dark Energy 
Experiment (HETDEX)^, and the WiggleZ surveys^ will 
measure the distance to higher redshifts within a few 
percent accuracy. 

^ http:/ /www. sdss3.org/cosmology.php 
^ http://hetdex.org/ 

^ http:/ /wigglez. swin.edu.au/Welcome. html 



The BAO signature appears as a small wiggle pattern 
in the galaxy power spectrum. Since the amplitude of 
BAO wiggle is very small a few percent), rather accu- 
rate theoretical models are needed. Especially, in order 
to determine the distance within a percent accuracy for 
the planned or ongoing surveys, we need to be able to 
predict the acoustic scale with much higher accuracies 
(~ 0.1%). However, there are complicated astrophys- 
ical processes such as the non-linear gravitational evo- 
lution, scale-dependent bias of galaxies, redshift space 
distortion, and the effect of massive neutrino. Many au- 
thors tackled these problems using numerical simulations 
(Meiksin et al. 1999; Seo & Eisenstein 2005; Huff et al. 
2007; Smith, Scoccimarro & Sheth 2007, 2008; Angulo 
et al. 2008; Takahashi et al. 2008; Seo et al. 2008; 
Nishimichi et al. 2009; Kim et al. 2009; Heitmann et 
al. 2009) and analytical perturbation theories (Crocce & 
Scoccimarro 2006, 2008; Jeong & Komatsu 2006, 2009; 
Nishimichi et al. 2007; McDonald 2007; Matarrese & 
Pietroni 2007, 2008; Eisenstein et al. 2007; Pietroni 
2008; Matsubara 2008a,b; Taruya & Hiramatsu 2008; 
Takahashi 2008; Nomura et al. 2008; Rassat et al. 2008; 
Sanchez et al. 2008; Padmanabhan & White 2008; Saito, 
Takada & Taruya 2009; Shoji, Jeong & Komatsu 2009; 
Taruya et al. 2009; Montesano et al. 2010). 

It is crucial to use not only accurate power spectra but 
also accurate covariance matrices in order to determine 
cosmological parameters from the galaxy power spec- 
trum. If the matter density fluctuations obey a Gaus- 
sian distribution, the covariance matrix has only diago- 
nal element and the relative error is simply given by the 
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square root of the number of modes in the survey area 
(e.g. Feldman, Kaiser & Peacock 1994). However, when 
the density fluctuations grow to the non-hnear regime, 
the mode couphng between different wavenumbers gen- 
erates non-zero off-diagonal elements, and the so-called 
non-Gaussian error is induced (e.g. Scoccimarro, Zal- 
darriaga & Hui 1999; Meiksin & White 1999). Rimes & 
Hamilton (2005,2006) first pointed out that there is little 
information contained in the power spectrum at quasi- 
nonlinear regime {k = 0.2 — 0.8/i/Mpc) due to the non- 
Gaussian error (see also Hamilton, Rimes & Scoccimarro 
2006; Neyrinck, Szapudi & Rimes 2006; Neyrinck & Sza- 
pudi 2007,2008; Lee & Pen 2008; Neyrinck, Szapudi & 
Szalay 2009; Lu, Pen & Dore 2009; Sato et aL 2009). For 
weak lensing (cosmic shear) analysis, the non-Gaussian 
error contribute the total error substantially (Cooray & 
Hu 2001; Sefusatti et al. 2006; Dore, Lu & Pen 2009; 
Takada & Jain 2009; Pielorz et al. 2009), but also may 
systematically shift the best fitting parameter (Hartlap 
et al. 2009; Ichiki et al. 2009). Especially, if there are a 
small number of parameters to be determined, the non- 
Gaussianity affects the errors significantly (see the dis- 
cussion in Takada & Jain 2009). 

In our previous paper (Takahashi et al. 2009, here- 
after T09), we used 5000 cosmological simulations to ob- 
tain the accurate covariance matrix of the matter power 
spectrum. This is a largest number of realizations ever 
done for the cosmological N-body simulation. We stud- 
ied the non-Gaussian error contribution to the signal-to- 
noise ratio for the measurement of the power spectrum, 
and found that the non-Gaussian error is important at 
small length-scale k > 0.2/i/Mpc. In this paper, we fur- 
ther investigate the non-Gaussian error contribution to 
the cosmological parameter estimation and the best-fit 
values using the likelihood analysis. We calculate the 
distribution of the best-fit parameters among the 5000 
realizations, and compare it with the results using the 
Fisher matrix analysis. We also study the distribution of 
the acoustic scale positions among the realizations. Our 
results in this paper can be used not only for the BAO 
analysis but also for the more general issue in the likeli- 
hood analysis of the non-linear matter power spectrum. 

Throughout the present paper, we adopt the standard 
ACDM model with matter density n„i — 0.238, baryon 
density fib — 0.041, dark energy density fi^ — 0.762 with 
equation of state w = — 1, spectral index Ug = 0.958, 
amplitude of fluctuations ag — 0.76, and expansion rate 
at the present time Hq = 73.2km s~^ Mpc~^, consistent 
with the 3-year WMAP resuhs (Spergel et al. 2007). 

2. MATTER POWER SPECTRUM AND ITS COVARIANCE 
MATRIX FROM NUMERICAL SIMULATIONS 

We follow the gravitational evolution of 256^ collision- 
less dark matter particles in a volume of 1000/i~"'^Mpc on 
a side using the cosmological simulation code Gadget- 
2 (Springel, Yoshida & White 2001; Springel 2005). 
We generate initial conditions following the standard 
Zel'dovich approximation using the matter transfer func- 
tion calculated by CAMB (Code for Anisotropics in the 
Microwave Background: Lewis, Challinor & Lasenby 
2000; also see Seljak & Zaldarriaga 1996). The initial red- 
shift is set to be z = 20. We use outputs at z = 3, 1 and 
0. To calculate the density fluctuations, we assign the 
N-body particles onto a 512^ rectangular grid using the 



cloud-in-cell scheme. Then we perform the Fourier trans- 
form and calculate the power spectrum in both real space 
and redshift space. We run 5000 Particle- Mesh(PM) sim- 
ulations to follow the non-linear evolution of the power 
spectrum and its covariance matrix in detail. We have 
checked that the power spectra of our simulations agree 
well with the result of the higher resolution TreePM sim- 
ulation, within 1(3)% for k < 0.2(0.4)ft./Mpc (here the 
Nyquist wavenumber is fc = 0.8/i/Mpc). If the initial 
redshift is set to be higher, z — 50, the results agree 
within 2(10)% for fc < 0.2(0 A)h/Mpc^. This is a suf- 
flcient accuracy for our purpose, which is to investigate 
the non-linear evolution of the power spectrum at BAO 
scales. 

Denoting Pi{k) as the power spectrum computed from 
the i-th realization, the ensemble averaged power spec- 
trum is estimated from the mean of the power spectra 
between 5000 realizations: 
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where Nr — 5000, the number of our realizations. Sim- 
ilarly, the covariance matrix between the spectra of fci 
and fc2 is estimated as 



cov(fci, fc2) = 



1 



^ i=i 



[/^.(fcl)-P(fcl)] [A(fe)-P(fc2) 



(2) 

The accuracy of the covariance is analytically estimated 
for the Gaussian density fluctuations (see Appendix). 
For example, the relative errors in the diagonal covari- 
ance terms are found to scale with the number of real- 
izations as our 5000 simulations provide a few 
percent accuracy. Clearly, our study achieves an unprece- 
dented accuracy of the covariance matrix estimation on 
BAO scales. 

The power spectrum covariance is formally expressed 
as a sum of the two contributions, the Gaussian and non- 
Gaussian terms (e.g. Scoccimarro, Zaldarriaga & Hui 
1999; Meiksin & White 1999): 

cov(fci,fc2) = ((P(fci)-P(fci)) (P(fc2)-P(fc2))) 

X r(kl,-kl,k'2,-k'2).(3) 

The first term arises from the Gaussian fiuctuations, 
while the second term is the non- Gaussian error arising 
from the mode coupling during the non-linear evolution. 
Here, P{k) — {P{k)) is the mean power spectrum, T is 
the trispectrum, the integral is done over the shell of the 
radius fci^2 with the width Afc in Fourier space, and Vk^ ^ 
is the volume of the shell given by Vfc = Ank^Ak. The 
expression in Eq.Q depends on the bin width, the first 
term is proportional to l/(l^Afc), while the second term 
is (X l/V. Hence, for the finer bin width, the impact of 
the Gaussian term becomes relatively enhanced. Note 

The agreement is achieved in real space. In redshift space, PM 
simulations somewhat underestimate the power spectrum at small 
scales (fc = 0.4/i/Mpc) by 20%. 
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however that the parameter estimation shown in the fol- 
lowing is independent of the bin width. Throughout this 
paper, the bin width is set to A/c = 0.01/i/Mpc.^ 

In this paper, we do not consider another non-Gaussian 
term in Eq. ([3]) arising from the finite survey volume (the 
so-called beat-coupling effect; Rimes & Hamilton 2006). 
Fluctuations with wavelength larger than the survey re- 
gion may contribute to the covariance on smaller scale. 
Although this effect can increase the covariance by over 
ten percent (T09), the main conclusions we draw in the 
present paper remain robust to the uncertainty. 

3. EFFECTS OF NON-GAUSSIAN ERROR ON LIKELIHOOD 
ANALYSIS OF COSMOLOGICAL PARAMETERS 

3.1. Parameter Estimation for Cosmological Parameters 

In this section, we study the effects of the non-Gaussian 
errors on the cosmological parameter estimation given 
the power spectrum measured from a hypothetical sur- 
vey of (1/i^^Gpc)'^ volume. We use the Fisher matrix 
formalism to estimate the accuracy of parameter estima- 
tion^ (Tegmark, Taylor & Heavens 1997): 
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aP(fci;x) 



dlnxi 



dP{k2 
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(4) 

where Xi denotes cosmological parameters, and the par- 
tial derivative such as dP/dxi is evaluated around the 
fiducial model. We include 5 parameters (therefore 
i = 1, 2, . . . , 5): the primordial power spectrum parame- 
ters,^ the normalization parameter Ag (not erg) and the 
spectral index rig, the baryon density fthh^, the dark mat- 
ter density i^ch"^, and the dark energy equation of state 
parameter w. We assume a fiat universe throughout the 
present paper. In Eq. Q we use Inxi (not Xi) as the 
variables such that the Fisher matrix gives the relative 
accuracy of a given parameter estimation: the marginal- 
ized error is then given as Axi/xi = {F^^)l/^ . Note that 
for ui, which has a negative value for the fiducial value, 
we simply compute wd In P/dw for the derivative. From 
Eqs.Q and Q, the estimation error Axi/xi is inversely 
proportional to the survey volume as Axi/xi oc V~^^'^. 

We need to compute the derivatives of the power spec- 
trum to compute the Fisher matrix in Eq. (j4]). For each 
cosmological parameter, we ran simulations with one pa- 
rameter slightly varied, while fixing other parameters to 
the fiducial values. We then compute the derivatives by 
the two-side differences of steps Axi/xi = ±0.05. We 
use 40 realizations are used for each parameter varia- 
tion. FigU] shows the derivatives of the power spectrum 
with respect to each cosmological parameter in real space 
(left column) and in redshift space (right column), re- 
spectively. The symbols are our simulation results: the 
red circles are for z = 0, the blue triangles for z = 1, 
and the green crosses for z = 3, respectively. The sensi- 
tivity of the power spectrum to cosmological parameters 

5 We also try the half bin-width, Afc = 0.005/i/Mpc, but our 
results in the following sections are almost same. 

® Here, we do not consider the Alcock-Paczynski effect (Alcock 
& Paczynski 1979) which would affect the measurement accuracy 
for the dark energy. 

We assume the primordial power spectrum given as Po{k) oc 
Aglk/ko)"" , where the pivot wavenumber ko is sot to ko = 
0.002/Mpc as employed by Komatsu et al. (2009). 



appears differently between in real- and redshift-space, 
due to the redshift distortion effects. For example, for 
a model with higher power spectrum normalization i.e. 
larger Ag, the power spectrum amplitudes are increas- 
ingly enhanced at larger k due to the stronger nonlin- 
earities in real space. In redshift space, however, the 
enhancement is significantly suppressed by the stronger 
finger-of-God effect, which arises from random velocity 
dispersion of dark matter particles in nonlinear objects. 
The relative amplitude of baryon acoustic oscillations is 
enhanced with increasing the baryon density. Finally, a 
change in w affects the power spectrum via the effect on 
the growth rate. We naively expect that the dependence 
of P{k) on w becomes scale-dependent in the nonlinear 
regime. However, the induced scale-dependence is weak 
and the dark energy parameter is very likely degener- 
ated with the galaxy bias in the measured galaxy power 
spectrum. 

The simulation results are compared with the analyt- 
ical predictions computed from the linear theory (the 
dashed curve) and the halofit (the short-dashed curve: 
Smith et al. 2003), respectively.^ The redshift-space 
power spectrum derivatives are compared with the linear 
theory. In the linear power spectrum in redshift space, 
we use the Kaiser approximation and do not include the 
Finger-of-God term. All the analytical predictions agree 
well with the simulation results at small k. In particu- 
lar, the halofit agrees with the simulations to within 10% 
accuracy. 

In our previous paper (T09), we studied the impact 
of the non-Gaussian error on the signal-to-noise ratio of 
power spectrum measurement: 
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,2<fem 



cov-i(fci,fc2)P(fci)P(fc2). (5) 



It was found that, in the linear regime, the S/N keeps in- 
creasing with increasing the maximum wavenumber fcmax 

as (S/N) oc fcmax- Howcver the S/N saturates at some 
fcmax in the weakly nonlinear regime, and stays nearly 
constant at larger fcmax > 0.2/i/Mpc due to the non- 
Gaussian errors. From these results one may naively 
guess that the parameter estimation is also significantly 
affected by the non-Gaussian errors when the power spec- 
trum information to the larger fcmax is included. In the 
following, we will study the impact of the non-Gaussian 
errors on the parameter estimation. 

Fig|2] shows the marginalized error on each parame- 
ter, Axi/xi = iF~^)l/'^, as a function of fcmax, where 
the power spectrum information up to a given fcmax is 
included. In each panel the symbols show the simulation 
results including the full covariance matrix. The simu- 
lation results are almost indistinguishable from the solid 
curves that are computed only by including the Gaussian 
error covariances, computed from simulations, in Eq. 
The agreement indicates that the Gaussian error assump- 
tion actually provides a good approximation for the pa- 

* We also compare the Lagrangian perturbation theory (LPT: 
Matsubara 2008a) with the simulation results in the manuscript 
in previous version (Takahashi et al. arXiv:0912.1381 vl). For 
the redshift-space spectrum, the LPT agrees well with the simula- 
tions at small k, but deviates significantly in the weakly nonlinear 
regime due to a too significant exponential damping, Pi^pTik) oc 
exp[— const. X k^] (Matsubara 2008a). 



4 



Takahashi et al. 



real space redshift space 




k (h/Mpc) k (h/Mpc) 

Fig. 1. — The derivatives of non-linear power spectrum with respect to cosmological parameters {As, Us, Qi,h^, Q,ch? , w) in real-space (left 
panel) and redshift-space (right). The cross, triangle and circle symbols show the simulation results at redshifts z = 1 and 0, respectively. 
The dashed curves are the linear theory predictions, the dotted curves are for the halofit. 



rameter estimation over scales of interest, even though 
the non-Gaussian errors have a significant impact on the 
S/N at fcmax > 0.2/i/Mpc. A more quantitative inter- 
pretation of these results will be given later. For com- 
parison, the dashed and dotted curves show the results 
obtained by using the linear theory and halofit to esti- 
mate the power spectrum as well as the Gaussian error 
covariances. These analytical predictions are far from the 
simulation results due to their inaccuracies in comparison 
with the simulations. Although the power spectrum and 
its derivatives in the halo fit agree within ~ 10% with 
the simulations (see FigH]), the parameter estimations 



are largely different as shown in FiglJl This is because 
when calculating the inverse matrix of the Fisher matrix, 
even a small errors in the Fisher matrix generates large 
errors in the inverse matrix. 

FiglHl shows the relative accuracies of each parame- 
ter estimation as a function of fcmax- There, we com- 
pare the results derived from the covariances with and 
without the non-Gaussian error contributions. The solid 
curves are the results where all the five parameters are 
included in the Fisher analysis, while the dashed curve 
shows the unmarginalized error on each parameter, i.e., 
the error is obtained by considering only one free pa- 
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rameter, Axi/xi = Fj; . In other words, the dashed 
curves correspond to the case where other parameters 
are well constrained by external data sets. The differ- 
ence between the solid and dashed curves is caused by the 
parameter degeneracies; the marginalized error becomes 
same as the unmarginalized error when the parameters 
are independent in the measured power spectrum. It is 
clear that, for the unmarginalized errors, including the 
non- Gaussian covariances degrades the parameter errors 
by a factor 4-5 for the redshift z = 0, and by a factor 
2-3 for z = 1, respectively. The level of the degradation 
is similar to that of the S/N as found in T09. There- 
fore, the impact of non-Gaussian covariance errors is sig- 
nificantly mitigated by the parameter degeneracies (see 
also, Neyrinck & Szapudi 2007; Takada & Jain 2009). 
The degradation of S/N increases the full Fisher ellipsoid 
volume, and then individual parameters are not tightly 
constrained due to the parameter degeneracies in such 
a high-dimension parameter space. As shown in FiglJl 
the non-Gaussian effect becomes negligible for the multi- 
parameter fitting. We comment that this conclusion is 
independent of the volume of the survey. 

In the upper two panels for and the short dashed 
curves show the results for the two parameter fitting case 
{As, Us), which are very similar to the solid curves. In 
reality, parameters that describe galaxy bias need to be 
further included. We thus conclude that the impact of 
the non-Gaussian errors is less important than the pa- 
rameter degeneracies, and that the Gaussian covariances 
can provide a good approximation to obtain the statisti- 
cal uncertainty of given parameters. 

3.2. Distribution of Best-Fit Parameters 

Nonlinear structure formation causes non-Gaussian 
distributions of the power spectrum estimators at small 
scales, as studied in, e.g., T09, in detail. Here, utilizing 
our 5000 realizations, we quantify the distribution of pa- 
rameter estimation taking into account the non-Gaussian 
covariances and the marginalization over other parame- 
ters. To this end, we simply use the x^-fitting analysis 
given as 

X?(x) = cov-i(fci,fc2) [P(fci;x) - A(fci) 

fcl,2<fcmax 



[P,(fci)-P(A:i,Xfid)] 



aP(fc2;x) 



where x — {As,ns,ftbh^ ,^ch^ 



X P(fc2;x)-P,(fc2)J ,(6) 

, Pi{k) is the power 
spectrum estimator of th i-th realization and P(fc; x) is 
its mean. For this analysis we simply use the real-space 
power spectrum. 

The variation in the power spectrum around the fidu- 
cial model can be expressed as 



P(fc;x) ~P(fc;Xfid) 



dP{k; x) 



dy 



(X - Xfid) 



(7) 



fid 



Recall that the best-fit parameters are estimated by min- 
imizing the ■ The best-fit parameters for the i-th real- 
ization can be estimated by inserting Eq. ([7]) into Eq. ([6]) 
(e.g. Huterer & Takada 2005; Joachimi & Schneider 
2009): 



(Xbt - Xfid), = 



F- 



E 



cov ^{ki,k2) 



.(8) 



fid 



Thus the best-fit parameters are generally different from 
the fiducial values depending on the distribution of Pi or 
how Pi deviates from the ensemble average expectation 
P at each wavenumber. Strictly speaking, we need to 
vary the covariance as a function of cosmological models 
in Eq. ([8]), but we here simply employ the covariance 
for the fiducial model assuming that the variations in 
the covariance are small (see Eifler, Schneider & Hartlap 
2009 for this issue). 

FigU shows how the best-fit values of w are dis- 
tributed among 5000 realizations. Note that, as can 
be seen in Eq. ([S]), the dark energy constraint includes 
the power spectrum amplitude information in addition 
to the BAO features. The left-panel shows the result 
for one parameter fitting {w), while the right panel for 
the three parameter fitting {w, As^ris), where the best- 
fit w is derived by including marginalization over the 
two parameters {AstUs). The latter corresponds to the 
case that the other parameters {VlJt? ,Q,bh?) are well 
constrained by external information such as the CMB 
and/or the Big Bang Nucleosynthesis (BBN). The red 
(black) curve shows the result obtained when the non- 
Gaussian errors in the Fisher matrix and the covariance 
in Eq. ([8]) are included (not included). The correspond- 
ing best-fit parameter deviations are plotted in the unit, 
[(w/wfid) — , where a is set to the Fisher errors with 
and without the non-Gaussian errors for the red and 
black curves, respectively. For example, the parameter 
deviations Kw/wfid) ~ ^ 3 mean that the parameter 
deviations are within ±3(7 confidence level regions. 

The distribution of the best-fit w looks nearly sym- 
metric: the nonlinear power spectrum does not shift the 
w-parameter to either of negative and positive sides from 
the fiducial value. The left panel (one-parameter fit- 
ting case) shows that including only the Gaussian er- 
rors makes the distribution of the best-fit w broader 
than the statistical confidence region. Clearly, a strong 
evidence on w ^ —1 may be incorrectly derived with 
high chances under the Gaussian assumptions. However, 
the red curves demonstrate that such apparent devia- 
tions can be corrected if we properly take into account 
the non-Gaussian errors for the statistical confidence re- 
gions. The right panel shows that, for a multi-parameter 
fit, the difference between the results with and with- 
out the non-Gaussian errors is significantly suppressed 
due to the parameter degeneracies. Note that our re- 
sults in FigH] is independent of the assumed survey vol- 
ume, because [w/{w)f^d — 1] oc V~^/'^ from Eq.® and 
a cx 

1^-1/2. Although we show the result only for w, 
essentially the same results are obtained for other pa- 
rameters {As^risT^bh^ T^ch^)- 

In FiglSl we quantify how the best-fit values of w are 
systematically different when including or ignoring the 
non-Gaussian errors. The horizontal axis is the differ- 
ence between the best-fit parameters, {wq — WNG)/wfid, 
divided by the Icr statistical confidence error derived by 
including the non-Gaussian errors. The left panel shows 
that the differences are smaller than the la confidence 
regions, and the right panels show even much smaller 
differences for the three-parameter fitting. Therefore we 
again conclude that the non-Gaussian errors do not cause 
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real space redshift space 




k^ax (h/Mpc) kn^ax (h/Mpc) 

Fig. 2. — The marginalized errors of each parameter when including the power spectrum information up to the maximum wavenumber 
fcmax in the horizontal axis. The left- and right panels show the results for the real- and redshift-space, respectively. The symbols in each 
panel are as for the previous figures. The solid curves which lie almost on top of the symbols show the results obtained only by including the 
Gaussian errors in the Fisher analysis. The agreement indicates that the Gaussian error assumption is a good approximation for parameter 
estimations even at small scales (see text for the details). The dashed and dotted curves are the analytical predictions that are derived 
using the linear theory and halofit for the power spectrum and the Gaussian covariance, respectively. 



any significant bias in the best-fit value compared to the 
statistical confidence regions including the non-Gaussian 
errors. 

4. ACOUSTIC PEAK POSITIONS 

A more robust method to constrain dark energy is us- 
ing the BAO peak positions. The BAO peak positions 
are characterized basically by one parameter, the stretch 
parameter (see below), and therefore the non-Gaussian 
errors may significantly affect the accuracy of the peak 



position determination. Here we use the distribution of 
the BAO peak positions obtained from our 5000 real- 
izations. Wc employ the method developed in Percival 
et al. (2007) to estimate the peak positions (see also 
Nishimichi et al. 2009). We first divide the measured 
P{k) by a smooth model, which is constructed by adopt- 
ing the cubic B-spline function to fit the binned power 
spectrum over a range of wavenumbers binned with the 
width Afc = O.Glft./Mpc. 
FiglHshows the power spectrum divided by the smooth 
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Fig. 3. — The ratio of the marginahzed errors with and without the non-Gaussian errors, as a function of fcmax. The solid curves show 
the results for our fiducial set of five cosmological parameters as shown in each panel. For comparison the dashed curves show the results 
for the unmarginalized errors or equivalently for the case of one parameter fitting in each panel. Also in the panels for As and ris the 
dotted curves show the results for the two parameter fitting of (As, n^), which appear to be between the solid and dashed curves. It is 
clear that the non-Gaussian errors degrade the unmarginalized errors by up to a factor 5. 



model: the data points are the average spectrum of 
5000 realizations and the errors the \a variation ranges, 
AP{k) = covi/2(fc,fc). Clearly, the BAG features are 
smoothed out at larger k and at lower redshifts due to 
stronger nonlinearities. 

To make parameter forecasts, let us define the ratio 
power spectrum, i?(fc), as 



R{k) 



P{k) 



smooth 



(fc) 



(9) 



Then we can introduce the stretch parameter a which 
characterizes a shift of the BAG peak phases via the 



transform k ak in R[k). The power spectrum with 
the stretch parameter a is given as 



P{k;a) = Psmootii(A:)i?(afc), 



(10) 



and a = 1 is the fiducial model. The Fisher information 
matrix for the stretch parameter a is computed as 

T? Sr -ifi. 1. ^ dP{ki;a) dP{k2]a) , . 

Faa^ COV (fci,fc2) — — . (11) 



da 



da 



Since we focus on the BAG peak locations, we treat only 
a as a free parameter, and hence the precision of deter- 
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Fig. 4. — The distribution of the best-fit parameter w estimated for the power spectra of 5000 realizations in real space. The red (black) 
curve shows the result obtained with (without) the non-Gaussian error covariance in Eq. l(8]l. The dotted, solid and dashed curves are for 
for femax = 0.1, 0.2 and 0.4/i/Mpc, respectively. The left panel shows the best-fit values of to for one parameter fitting [w alone), while the 
right panel shows the results for three parameter fitting (to, j4s,ns). The horizontal axis is defined as (to/tOflj — l)/^, where a denotes the 
1(T confidence regions computed from the Fisher matrix with and without the non-Gaussian errors for the red and black curves, respectively. 



parameter 




Fig. 5. — As for the previous figures, but for the distribution of the differences between the best-fit values of w computed from Eq. Jsjl 
with and without non-Gaussian errors. The (Tnq is the \a error computed from the Fisher matrix with the non-Gaussian errors. 



mining a for the given power spectrum measurement is 
given as Aa = baa ■ 

FiglT] shows the Icr error, Aa, as a function of fcmax 
up to which the power spectrum information is included. 
The symbols are the results including the non-Gaussian 
errors in the Fisher analysis, while the solid curves are 
for the Gaussian errors. The Gaussian error assump- 
tion appears to be valid even for large A:niax- This is 
because the BAO features are erased at the weakly non- 
linear scales fc > 0.2ft,/Mpc, where the non-Gaussian 
errors are more significant. There is little information 
on the acoustic scale at the nonlinear scales. At 2; = 
the accuracy improves significantly around the first peak 
(fc 0.06/i/Mpc) and the second peak (fc ~ 0.12/i/Mpc). 
Hence almost all the information on the acoustic peaks 
are obtained for k < 0.15/i/Mpc at z = 0. From Figd 
Aa ~ 1% can be achievable for a survey with (Gpc//i)'^ 



volume coverage. 

Finally, we investigate the distribution of best-fit a 
among 5000 realizations. Given the power spectrum 
measurement for the «-th realization, the for estimat- 
ing a is given as 

X?(") = X! COV~^(fci, fc2)Psmooth(fcl)^'smooth(fc2) 



2<fcn 



(12) 



Here, we have used the power spectrum Pi{k) for the i- 
th realization, not the mean P{k), to obtain the smooth 
power spectrum in the denominator of the ratio Ri(k), 

Fig IS] shows the distribution of the best-fit shift pa- 
rameters a for the cases of fcmax = 0.2 and 0.4/i/Mpc. 
The horizontal axis is a — 1 divided by the la error ctq, 
obtained by the Fisher matrix. The red curves are the re- 
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Fig. 6. — The real-space power spectrum divided by the 
smoothed spectrum for 2; = 0, 1 and 3, respectively. The symbols 
denote the mean of the power spectra among 5000 realizations of 
(1/i/Gpc)^ volume, while the error bars denote the la scatters. 
The solid curve is for the linear P{k). The BAO features are more 
erased at higher k and at lower redshifts due to stronger nonlin- 
earities. 



Fig. 8. — The distribution of the shift parameter a for fcmax = 0.2 
(solid curves) and 0.4/i/Mpc (dashed curves), respectively. The 
horizontal axis o — 1 is divided by the Icr error computed from the 
Fisher matrix. The red and black curves are the results with and 
without the non-Gaussian errors, respectively. The distribution of 
a is well described by a Gaussian distribution with the width of 
the la Fisher error. 
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Fig. 7. — The la error on the BAO peak position parametrized by 
the stretch parameter a (see text for the definition). The symbols 
are the results computed from simulations, while the solid curves, 
which are almost top on the symbols, show the errors computed 
assuming the Gaussian errors. 
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Fig. 9. — The distribution of the differences between the best-fit 
values of a with or without including the non-Gaussian covariance. 
The solid (dashed) curve is for fcmax = 0.2(0.4)/i/Mpc. 



suits including the non-Gaussian errors, while the black 
curves are for the Gaussian errors. The two results are in- 
deed very similar. The distribution is well described by a 
Gaussian function with the width given by the Icr Fisher 
error. A recent work by Seo et al. (2009) compares the 
distribution of the acoustic scales with the result from 
the Fisher matrix analysis. Although their analysis is 
slightly different from ours (they use a fitting formula for 
the non-linear P{k) when calculating the acoustic scale), 
their results are broadly consistent with the results shown 
here. 

Fig|9]shows the difference between the best-fit shift pa- 



rameters when including or ignoring the non-Gaussian 
errors in Eq. (fT2|) . The difference is smaller than ~ 0.2a. 
Hence we conclude that the non-Gaussian covariance 
does not cause a substantial systematic error in the BAO 
peak determination. 

5. DISCUSSION AND CONCLUSION 

We have studied the effects of the non-Gaussian er- 
ror on the parameter estimations and the distribution 
of the best-fit parameters for the cosmological param- 
eters (in section 3) and the acoustic scale position (in 
section 4). We have found that the non-Gaussian error 
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is important for the parameter errors if there is only one 
fitting parameter. The measurement error degrades up 
to factor 5 for fcmax = 0.4. However, if there are more 
than two parameters, the impact of the non-Gaussian 
errors are insignificant due to severe parameter degen- 
eracies in the matter power spectrum. For the acoustic 
scale, the non-Gaussian errors do not affect the acoustic 
scale determination, even though there is only one fit- 
ting parameter. This is because that the acoustic scale 
is determined mainly by the linear scales where the non- 
Gaussian covariance is not important. 

Throughout this paper, we discussed the covariance 
matrix of the matter power spectrum. However, for real 
galaxy survey, we should include the effects of the halo 
and galaxy bias. However addressing this issue is not 
easy at present, because we need the large number of real- 
izations of high resolution simulations including halo and 
galaxy formation. For large scale {k < 0.1/i/Mpc), the 
covariance of the halo power spectrum is consistent with 
the Gaussian error with the shot noise term (Angulo et al. 
2008; Smith 2009). Because, in the linear regime, only 
the sample variance dominates the covariance. For small 
scale {k > 0.1/i/Mpc), Smith (2009) recently showed that 



the covariance is larger than the Gaussian error predic- 
tion and the higher mass halo have stronger covariance 
due to the non-linear gravitational mode coupling. He 
compared the cluster-sized halos in the two mass ranges 
(M > 10"Mo and 1 x IO^^Mq < M < 2 x IO^Mq), 
and hence we expect the galactic halo (M < 10^^ Mq) 
would have a weaker covariance. 

In our previous paper (T09), we compare the power 
spectrum covariance in numerical simulations with an 
analytical models such as perturbation theory and halo 
model. We calculated the diagonal and off-diagonal 
terms of the covariance matrix and the signal-to-noise 
ratio in the both models, and found that the halo model 
reasonably well reproduces the simulation results. Sev- 
eral authors also compare them and reached the same 
conclusion (e.g. Cooray & Hu 2001; Neyrinck et al. 2006; 
Neyrinck & Szapudi 2007; Sato et al. 2009). 

Our simulation results of the 5000 power spectra Pi{k), 
the derivative of P{k) with respect to the cosmologi- 
cal parameters, dP{k)/dxj, and the covariance matrix 
cov(fci, k2) are available as numeric tables upon request 
(contact takahasi@cc.hirosaki-u.ac.jp) . 



APPENDIX 

VARIANCE OF THE COVARIANCE MATRIX FOR THE GAUSSIAN DENSITY FLUCTUATIONS 
The covariance matrix estimated from realizations is given by, 

cov(fci, k2; Nr) = ^Yl [^i(^i) - ^(^i)] [-^^(^2) - P{k2)] , 

^ i 

= l^E^^(W(/^2)-^i:A(fcl)l^E^:'(^2). 



The variance of the covariance is given by, 

var[cov(A;i,A;2;iVr)] = (cov2(fci, fe; iVr)) - (cov(fci, ^2; iV^))^ 

= ^HliPiiki)Pi{k2)Pj{ki)Pjik2)) - {Pi{k,)h{k2)){m,)p^{k2)) 

[(A(fcl)A(fe)P,(fcl)Pfc(A2)) - (A(fcl)A(fc2))(P,(ftl)Pfc(fc2)) 
+ E [{Piikl)Pjik2)Pk{kl)Pl{k2)) - {Pi{kl)Pj{k2)){Pk{k,)Pl{k2)) 



(Al) 



(A2) 



i,j,k,l 



The above equation further reduces to 





1 








2 - 








1 


+ 








2 


+ 







(A3) 



{P\ki)P{k2)){P{k2)) - {P{k^)P{k2)){P{ki)){P{k2)) + (fcl ^ k2) 
[{P\k,)) - {P{k,)f) {P{k2)f + (fcl ^ fc2)' 
(P(fcl)P(fc2)) - (P(fcl))(P(fc2))] (P(fcl))(P(fc2)) 

Here we ignored the terms of the order of ^'^d higher. 

For the Gaussian density fiuctuations, the n-th moments (P") can be obtained using the probability distribution 
function of P{k) (the chi-squared distribution function, see Eq.(Bl) in T09). For the diagonal parts, we have 



var[cov(fci,fci;7V,)] = — (^-^ + 



P\k) 



(A4) 
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Fig. 10. — The dispersions among the power spectrum covariances each of which is estimated from the A^r reahzations (a subset of the 
while 5000 realizations), as a function of Nr. In the vertical axis, the dispersion is defined in Eq. IIA7I I. The panel shows the result for the 
off-diagonal parts for varying the wavenumber bins and the bin widths. The color symbols are the simulation results, while the solid curves 
denote the theoretical prediction for the Gaussian density fluctuation. The plots explicitly show that the power spectrum covariances are 
estimated at a sub-percent level accuracy by using our whole 5000 realizations. 



Hence we have 

var [cov(fci,fci; iVr)] _ _2_ , . 

cov2(fci,fci) ~ ^ ' 

which is consistent with our numerical finding in our previous paper (see left panel of Fig. 12 in T09). Hence if we need 
10%(5%) accuracy in the covariance, we have to prepare 200(800) realizations. 
For the off-diagonal parts we have 

var[cov(fci,fc2;^r)] = ]J^]^]|;^'(^i)^'(^2) (A6) 

Let us define the relative errors as, 

^2 ^ var[cov(fci,fc2;^r)] ^-^^s^ 

'^"^ COv(fci, fci) C0v(fc2, fc2) 

In FigfTUl we show the relative errors a'^^^ in Ea. (|A7p as a function of N,-. The solid line is the theoretical prediction 
in Eg. (I ATI) , and the symbols are the results of our numerical simulation for various scales (fc = 0.05, 0.2, 0.4/i/Mpc) 
and bin-width (Afc — 0.005, 0.01/i/Mpc). Analytical results in Eq. (jA6p fit our simulation data well. 
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